Adaptation Reduces Variability of the Neuronal Population Code 
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Sequences of events in noise-driven excitable systems with slow variables often show serial correla- 
tions among their intervals of events. Here, we employ a master equation for generalized non-renewal 
processes to calculate the interval and count statistics of superimposed processes governed by a slow 
adaptation variable. For an ensemble of spike-frequency adapting neurons, this results in the reg- 
ularization of the population activity and an enhanced post-synaptic signal decoding. We confirm 
our theoretical results in a population of cortical neurons recorded in vivo. 

PACS numbers: 87.19.1s, 05.40.-a, 87.19.1o, 87.19.1c 



Statistical models of events assuming the renewal prop- 
erty, that the instantaneous probability for the occur- 
rence of an event depends uniquely on the time since 
the last event, enjoys a long history of interest and ap- 
plications in physics. However, many event processes in 
nature violate the renewal property. For instance, it is 
known that photon emission in multilevel quantum sys- 
tems constitutes a non-renewal process [l|. Likewise, the 
time series of earthquakes typically exhibits a memory of 
previous shocks [2j, as do the times of activated escape 
from a metastable state, as encountered in various sci- 
entific fields such as chemical, biological, and solid state 
physics 3]. Often, the departure from the renewal prop- 
erty arises when the process under study is modulated 
by some slow variable, which results in serial correlations 
among the intervals between successive events. In par- 
ticular, the majority of spiking neurons in the nervous 
systems of different species show a serial dependence be- 
tween inter-event intervals (ISI) due to the fact that their 
spiking activity is modulated by an intrinsic slow vari- 
able of self-inhibition, a phenomenon known as spike- 
frequency adaptation |4|. 

In this letter, we present a non-renewal formalism 
based on a population density treatment that enables 
us to quantitatively study ensemble processes augmented 
with a slow noise variable. We formally derive general ex- 
pressions for the higher-order interval and count statistics 
of single and superimposed non-renewal processes for ar- 
4i bitrary observation times. In spiking neurons, intrinsic 

42 mechanisms of adaptation reduce output variability and 

43 facilitate population coding in neural ensembles. We con- 

44 firm our theoretical results in a set of experimental in vivo 

45 recordings and analyse their implications for the read-out 

46 properties of a postsynaptic neural decoder. 

47 Non-renewal Master Equation. We define the limiting 

48 probability density for an event given the state variable 

49 x by the so-called hazard function h x {x,t) where t de- 

50 notes explicit dependence on time due to external input 

51 following [3, [6( . Here, we assume x has a shot- noise-like 

52 dynamics, which is widely used as a model of spike in- 



duced neuronal adaptation 



<t)/T + q^2 k 6(t-t k ), 



(1) 



54 where 5 is the Dirac delta function, tk is the time of the 

55 k th event, and q is the quantile change in x at each event. 

56 The dynamics of x deviates from standard treatments of 

57 shot-noise (such as in [7J ) in that the rate of events has 

58 a dependence on x as expressed by the hazard function 

59 h x {x,t). It is straightforward to show that the distribu- 

60 tion of x in an ensemble, denoted by Pr(x, t), is governed 
6i by 



d t Pr(x, t) = d x [- Pr(i, t)] + h x (x - q, t) Pr(x - q, t) 



h x (x,t)Pi(x,t). 



(2) 



62 Much insight can be gained by applying the method of 

63 characteristics [8[ to establish a link between the state 

64 variable x and its time- like variable t x . For Eq. ([T]) 

65 we define t x = rj{x) := —rhi(x/q), whereby 4rt x = 1. 

66 When an event occurs, t x i-> ip(t x ), where ip(t x ) = 

67 ^{rf 1 ^) + q) = — rln(e~ ta; / T + 1) with its inverse 

68 given by vb{t x )~ l = — rln(e~* x / r — 1). Thus, we de- 

69 fine h(t x ,t) := h x (j] (t x ),t). This transformation of 

70 variables to t x elucidates the connection of the model to 

71 renewal theory. Here, the reset condition after each event 

72 is not t x t— > (renewal) but t x h- > r](x + q) [5|- Therefore, 

73 the variable t x that we may call a 'pseudo age' is a general 

74 state variable that no longer represents the time since the 

75 last event (age). Transforming variables in Eq. @ from 

76 x to t x yields in the steady state 

dt, Pr(0 = -h(t x ) Pr(t x ) 

+ (1 - 9 (£ x )) [ft^- 1 (4) Pr(V>- X (**))], (3) 

77 where @o(t x ) is the Heaviside step function, and for con- 

78 venience we defined V ;_1 (^ > 0) = 0. An efficient ai- 
re gorithm for solving Eq. ([3]) is given in [f|. We denote 
so this solution by Pi eq (t x ). Further, the time- like trans- 
si formation in Eq. <j3j> allows computation of the ISI by 

82 analogy to the renewal theory [6| and also permits the 

83 comparison to the master equation for a renewal pro- 

84 cess as given in Eq. (6.43) in [9(. The distribution 



85 of t x just prior to an event is a quantity of interest 

86 and it is derived as Pr*(t x ) = h(t x )Pr eq (t x )/r eq , where 

87 r eq = J h(t x )Pr eq (t x )dt x is a normalizing constant and 

88 also the process intensity or rate of the ensemble. Simi- 

89 larly, one can derive the distribution of t x just after the 

90 event, Pr 1 "^) = Pr*^ -1 ^))^ -1 ^ 0- Then the 

91 relationship between t x and the ordinary ISI distribution 

92 can be written as 



128 (3.3) in [ll|], the second moment of the count statistics 

129 can be derived. Thus, we obtain the Fano factor 

J T = 1 + (2/T) J T (T - u)A(u)du - r eq T, (8) 

no The asymptotic property of F = limx-s-oo Jt can be de- 
i3i rived from the result stated in Eq. (7.8) in 



/+oo 
h(t x + A)Q(t x + A)Pr f (t x )dt Xl 
-OC 



(4) 



93 where fl(t x + A) = e~ $t K** +«)<*". n ow the n th moment 

94 p n of the distribution and its coefficient of variation C v 

95 can be numerically determined. 

96 Counting Statistics. In order to derive the count dis- 

97 tribution, we generalize the elegant approach for deriving 

98 the moment generating function as introduced in |10| : let 

99 p n {tn, t x \t x ) be the joint probability density given its ini- 
ioo tial state t x , where t n stands for time to n th event and t " 
ioi is the corresponding adaptive state of the process. Thcre- 
102 after, one can recursively derive 

p n+1 {s,t x + i \ti)= [ p n ( S ,t x \t x )p(s,t: +i \t x )dt:, (5) 



103 where p n+1 (s,t™ +1 \t x ) = £[p n+1 (t n+1 ,t™ +1 \t° x )] and C is 

104 the Laplace transform with resepect to time, assuming 

105 pi(s,t x \t x ) = p(s,t x \t x ) [10j. Next, defining the opera- 

106 tor P„ (s) and applying Bra-Kat notation as suggested 

107 in [lCf, leads to the Laplace transform of n th events or- 
loa dinary density 

p n {s) = (1 | P„(s) | Pr f ) = (1 | [P(.s)]» | Pr f ), (6) 



where the operator P associated with p(s), which inter- 
estingly corresponds to the moment generating function 
of the sum of n non-independent intervals f n (s) as de- 
fined in 11|. Now, following Eqs. (2.15) in |llj Laplace 



transform of count distribution denoted as P(n, s). 

The Fano factor provides an index for the quantifica- 
tion of the count variability. It is defined as Jt = o~ T / pt, 
where a T and px are the variance and the mean of 
the number of events in a certain time window T. It 
follows from the additive property of the expectation 
that pt = L r(u)du and assuming constant firing rate 
Pt = r eq T. To calculate the second moment of P(n,s), 
we require A s = J2k Pk(s), thus 



4, = (l|P(«)(I-P(s))- l |Pr 1 



(7) 



where I is the identity operator. Note, assuming a re- 
newal interval distribution in Eq. (j4]) one obtains A r s = 



124 p(s)/(l - p(s)) and L 1 [r eq A s 



jA{u) is the joint 



density of an event at time t and another event at time 
t + u. Thus, the autocorrlation of events is A(u) — 
r eq [8{u) + A(u)]. Now, by using Eq. ([7]) and the Eq. 



11 



]im s _, [A s - l/(p lS )} = C t 2 [l/2 + Er=i ffc] - 1/2, (9) 

where ^ is the linear correlation coefficient between two 
k lagged intervals. Provided the limit exits, we find F = 

C?[l + 2E£i&]in|ia. 

Superposition. We now generalize our results on the 
counting statistics to the superposition of independent 
point processes. This is of practical interest in all cases 
where we observe superimposed events that stem from 
multiple independent process, e.g. in photon detection 
devices, or in the case of a postsynaptic neuron that re- 
ceives converging inputs from multiple lines. We study 
the superposition of k stationary, orderly, and indepen- 
dent processes. The ensemble process will have a rate 
f = Yll=i r i and following Eq. (4.18) in [HI A{u) = 
f + f^ 1 X)i=i r i[-A-i(u) —ri\. Here, for the sake of simplic- 
ity, we derive the desired relationship between C% and 
the ensemble F for k identical processes. To this end, 
we plug f and £[.A(u)] into the Eq. © and therefore 
it becomes lim s ^ [A - 1/(mi«)] = CV 2 [l/2 + S] - 1/2, 
where CV and H = X)i*li "« are t ne coefficient of vari- 
ation and the interval correlations of the superimposed 
process. Note that the left hand side of this equation and 
Eq. ^) arc simular. Thus, we obtian 



CV 2 fl 



2S] = C t 2 [l + 2E, C : i &]- 



(10) 



154 The left hand side of Eq. (fT0j) is indeed the Fano factor 

155 F of the ensemble process as desired. 

156 gests as k — > oo, CV" — > 1. 



Now, [13j sug- 
Interestingly, if all individual 

157 processes fullfill the renewal condition, it follows from 

158 Eq. (HO]) that F = C 2 = [1 + 25], and therefore if C 2 ^ 1 

159 the population activity is non-renewal with S < (S > 0) 
wo for processes with C 2 > 1 (C 2 < 1). This important find- 
i6i ing explains the numerical observation in [14J regarding 

162 emergance of non-renewal processes as the result of the 

163 superposition operation. 

164 Adaptation in a Neuronal Ensemble. In [6| it has 

165 been shown by an adiabatic elimination of fast vari- 

166 ables that the master equation description of a detailed 

167 neuron model including voltage dynamics, conductance- 
's based synapses, and spike-induced adaptation reduces to 

169 a stochastic point process simular to Eq. ([3]). The corre- 

170 sponding hazard function can be approximated as 



h x (x) = a t exp(-b t x), 



(11) 



171 where at and bt are determined by the time dependent 

172 statistics of inputs [5( and the equilibrium rate consis- 

173 tency equation r eq w h x (r eq qT) [6( with the solution 



W{abq T )/(bqT), 



112) 
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Experiment 




20 0.1 

Expected count in T 



FIG. 1. Adaptation reduces the Fano factor of the ensemble 
process. Left Magenta: Jt for arbitrary observation time T 
according to Eqs. ([8| and (|lip with bq = 1.4, a — 5.0 and r = 
400ms. Blue: Fano factor for equivalent renewal ensemble 
process with interval distribution of Eq. ((4]). Square Dots: 
Numerically estimated Fano factor for superposition of the 5 
realization runs of the full-detailed adaptive neuron model as 
in [y|. Dash-dotted line: C%- Right Magenta: Empirical Jt 
estimated from the pooled spike trains of 5 cortical neurons. 
Blue: Fano factor for the pool of shuffled spike trains. Dash- 
dotted: Average Cl of the 5 individual spike trains. 



209 noise in the neuronal population rate code. Our analysis 

210 of a set of cortical data suggests a reduction of > 50% 
2n for long observation times. The reduction of Jt in Fig. Q] 

212 becomes significant even for small observation times of 

213 s=a 2 average intervals, which is a relevant time scale for 

214 the transmission of a population rate signal. This result 

215 is reminiscent of an effect that has previously been ac- 

216 knowledged as noise shaping and weak stimuli detection 

217 expressed in the reduction of the low frequency power in 
2i8 a spectral analysis of spike trains with negative serial in- 

219 terval correlations [171 ] . Our result confirms their findings 

220 at the population level. 

221 Our second argument is concerned with the transmis- 

222 sion of a population rate signal. We may define a func- 

223 tional neural ensemble by the common postsynaptic tar- 

224 get neuron that receives the converging input of all en- 

225 semble members. To elucidate the postsynaptic effect 

226 of adaptation we simplify the ensemble autocorrelation 

227 function A(u) following [18j with an exponential approx- 

228 imation 



174 where W is the Lambert function. In the case of van- 

175 ishing adaptation (bq — > 0) the process apporaches the 

176 Poisson process with r eq — > a. 

177 We show in [5[ that the adaptation dynamics in Eq. ([T]) 

178 produces negative serial correlations £& < 0. The 

179 strength of serial correlation decays with increasing lag 
lso k and depends on the mean adaptation, E[x\ — r eq qr. 
i8i Such a vanishing of negative serial interval correlations 

182 with increasing lag is well supported by a large body of 

183 experimental evidence [J]. The departure from the re- 

184 newal property induced by adaptation reduces the Fano 

185 factor Eq. §£$ for the single process as well as for the 

186 population model of superimposed pocesses. 

is? We validate our theoretical result of the reduced Fano 

188 factor in a set of experimental spike trains of N = 5 

189 in vivo intracellular recorded neurons in the somatosen- 

190 sory cortex of the rat. The spontaneous activity of each 

191 of these neurons shows negative serial interval correla- 

192 tions [15[ where the empirical sum over correlation co- 

193 efficients amounts to an average J2 i=1 S,i = —0.28. We 

194 construct the population activity by superimposing all 

195 5 spike trains. Thereafter, we estimate the Fano factor 

196 as a function of the observation time and compare it to 

197 the case where, prior to superposition, renewal statistics 

198 is enforced for each individual neuron through interval 

199 shuffling. Our experimental observation in Fig. [1] (Right) 

200 confirms the theoretical prediction of a reduced Fano fac- 

201 tor simular to individual neurons [16j in the population 

202 level. 

203 Benefits for Neural Coding. We provide three argu- 

204 ments that demonstrate how the mechanism of spike- 

205 frequency adaptation benefits neural processing and pop- 

206 ulation coding. First, our result of a reduced Fano factor 

207 F < C^ for the population activity of stationary adap- 

208 tive processes (bq > 0) directly implies a reduction of the 



A{u) = r eq 5(u) + [(F - 1)/2t c ] exp(-u/r c ), (13) 

229 where the second term is the approximation of r eq A(u). 

230 For given observation time window u, and r c the reduc- 

231 tion of F implies that A r u < A u . Therefore, the postsy- 

232 naptic neuron receives inputs from an adaptive ensemble 

233 that expresses an extended autocorrelation structure as 

234 compared to the inputs from a non-adaptive ensemble. 

235 Following the theory on the effect of input autocorrcla- 

236 tion on signal transmission in spiking neurons as dcvel- 

237 oped in [18| . a longer r c reduces the input current fluc- 

238 tuations and this facilitates a faster and more reliable 

239 transmission of the modulated input rate signal by the 

240 postsynaptic target neuron. 

241 Finally we argue that a postsynaptic neuron can bet- 

242 ter decode a small change in its input if the presynaptic 

243 neurons are adaptive. To this end, we compute the in- 

244 formation gain of the postsynaptic activity, between two 

245 counting distributions of an adaptive presynaptic enscm- 

246 ble when h x {x) is adiabatically transferred to h x (x — e) 

247 with a small change e in the input ensemble. It is con- 

248 venient to use p„(s) which associated with counting dis- 

249 tribution P(n,s). Thus, we apply the Kullback-Leibler 

250 divergence to the Eq. © before and after the adiabatic 

251 change in the input 

D K L(p e n \\p n ) = EiPi(s)HPi(s)/Pi(s)). (14) 

252 Using Eq. we obtain D KL (p e n \\p n ) = A e s [\n(A e s /A s )}. 

253 Due to Eqs. |T]) and (fT2"]) . the mean adaptation after the 

254 change is E[x e ] = rqr\ . If e > it follows that r\ q > r eq . 

255 Therefore the mean adaptation level increases and the 

256 adapted process exhibits stronger negative serial corre- 

257 lations and A\ > A s . Thus, by Eq. (fT3|) . it is straight 

258 forward to deduce that Dkl > D r KL . for renewal and 

259 adaptive processes with identical interval distributions. 




r (Hz) 
eq x ' 



FIG. 2. Information gain per spike due to adaptation. Left: 
Transfer of equilibrium rate for fixed e change of the input 
in adaptive and Poisson model. Right: Kullback-Leibler Di- 
vergence per extra spike as the measure of information gain 



300 ter equation [6| assumes incoherent input fluctuation as 

301 in the mean-field theory, where common input is negli- 

302 gible. Treating a network with coherent fluctuations as 

303 encountered in finite size networks requires an alternative 

304 derivation of the hazard function [5[ . 
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th 



for n event density of adaptive and Poisson processes while 



200ms and e = O.OlnS with the same initial 



and 



260 We now compute the information gain of the adap- 
26i tivc ensemble process relative to a matched Poisson rate 

262 model. For different initial rate values r eq we assume a 

263 small but fixed increase e in the input that we express in 

264 parameter changes a e and b e in Eq. (jlip as outlined in 

265 |5|. This leads to an increase k = r c — r eq in rate that 

266 is effectively constant over a wide range of initial values 

267 r eq (Fig. 2, Left). In the rate model, assuming the same 

268 initial value of r eq , the same input step leads to a higher 

269 equilibrium rate increase K Pmsson > Ki which depends 

270 on the inital rate (Fig. 2, Left) because the rate model 

271 lacks a mechanism of self-inhibition, which in the adap- 

272 tive model counteracts the rate increase. Thereafter, we 

273 compute the Kullback-Leibler divergence for both mod- 

274 els and normalize it by the change in the output rate k. 

275 The result in Fig. 2 (Right) shows that Dkl/k is larger 

276 for the adaptive model than for the rate model across 

277 the range of tested input rates. Thus, the information 

278 per extra spike is larger in the adaptive ensemble than in 

279 the renewal ensemble, and a postsynaptic neuron can dis- 

280 criminate small changes e more efficiently, even though 
28i the absolute change in firing rate is lower. 

282 Discussion. Our results point out a new aspect of spike 

283 frequency adaptation that benefits the reliable transmis- 

284 sion and postsynaptic decoding of the neural population 

285 code. This aspect adds to the known properties of com- 

286 pression and temporal filtering of sensory input signals 

287 [19( in spike frequency adapting neurons. The specific 

288 result of Eq. (fTU)) is also of practical consequence for the 

289 empirical analysis of the so-called multi-unit activity. By 

290 estimating Fano factor and serial correlations we readily 

291 obtain an estimate of the average C v and serial correla- 

292 tion of the individual processes. 

293 We developed a new formalism to treat event emit- 

294 ting systems that arc influenced by a slow state variable, 

295 and we provide a number of useful general results on 

296 the higher order event statistics of superimposed renewal 

297 and non-renewal event processes, which are applicable 

298 to a wide range of event- based systems in nature [fj. 

299 The derivation of the state dependent hazard and mas- 
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